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1. INTRODUCTION 


In spectral methods, the solution of a partial differential equation is represented by a 
truncated expansion of eigenfunctions of a singular Sturm-Liouville problem. This choice is 
responsible for the superior approximation properties of spectral methods. This is, of course, 
most evident for problems possessing smooth solutions in which case the series expansion 
converges faster than any inverse power of n as n — » oo. This phenomenon is known as 
spectral convergence. The method of calculating the expansion coefficients determines the 
type of spectral approximation: Galerkin, tau, or collocation. Only collocation methods are 
considered in this paper since they are applicable to a wide class of problems. 

Spectral methods are most easily applied to problems defined in rectangular or circular 
regions in which case Chebyshev or Fourier series, respectively, are appropriate. However, 
the natural choices of expansion functions for a problem defined in an irregular geometry 
are unwieldy and inefficient to use and need to be computed for each new irregular region 
[12]. There are two ways of overcoming these difficulties, namely, mapping and patching [2, 

4 . 12 ]. 

The mapping technique involves transforming the irregular region into a simpler one 
by using a coordinate transformation. Spectral methods can then be applied in the simpler 
region using standard expansion techniques [12, 16]. The patching method divides the region 
into a number of simpler subregions or elements. A spectral approximation to the solution of 
the differential equation is sought within each element. The representations are patched by 
imposing continuity conditions across interfaces. This results in a coupling of the expansion 
coefficients in contiguous elements. Spectral domain decomposition methods combines the 
flexibility of the finite element method with the superior approximation properties of the 
spectral method [3, 13, 14]. 

A number of different spectral domain decomposition techniques have appeared in the 
literature [3, 5, 7, 10, 13, 14]. The main differences between these variants lie in the choice 
of trial functions and the treatment of the continuity conditions at element interfaces. In 
the spectral element method [13] conforming elements are used and C l continuity across the 
interfaces is achieved implicitly through a variational principle. The global element method 
[3] uses trial functions which are nonconforming. A modified functional is used to ensure that 
the interface continuity conditions are satisfied. In the present paper, we advocate the use 
of conforming spectral domain decomposition techniques and describe a collocation strategy 
for achieving this. 

As our test problem, we consider Stokes flow through an abruptly contracting channel 
with contraction ratio 1 : a. A conforming spectral domain decomposition of this geometry 
divides the flow region into three rectangular semi-infinite subdomains with common point 
(0, a) as shown in Figure 1. In previous work [5, 15], the authors consider nonconforming sub- 
regions because of their ease of implementation (see Figure 2). Although this strategy works 
well for the Stokes problem [5], a lack of interface continuity appears for the Navier- Stokes 
problem [6] at moderate values of the Reynolds number, eventually causing the method to 
break down. The resulting spectral approximations are not pointwise continuous across the 
subregion interfaces. In this paper, we alleviate any possibility of discontinuous solutions and 
normal derivatives across the interfaces by using a carefully constructed collocation scheme. 

Efficient direct methods for solving the collocation equations in rectangularly decompos- 
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able domains are described in [15] for the two subdomain example. This method, which is 
basically a block Gauss method, is applied here to efficiently invert the coefficient matrix. 
The matrices resulting from a spectral collocation discretization assume a block structure 
in which the non-zero blocks are full. Economically viable solution techniques need to take 
advantage of this matrix structure. 

In Section 2, we define the Stokes problem in the channel contraction and derive the 
boundary conditions in terms of the stream function. In Section 3, we describe a domain 
decomposition of the flow region and the spectral approximation within each subregion. 
The collocation strategy which gives pointwise C 1 continuity of the stream function across 
the subregion interfaces is developed in Section 4. The solution of the spectral collocation 
equations using the capacitance matrix technique is described in Section 5. Numerical results 
are presented in Section 6 and concluding remarks are made in Section 7. 


2. The Governing Equations 


The governing equations for the planar inertialess flow of an incompressible Newtonian 
fluid assume the mathematical form 

V ■ v = 0, (2.1) 

V • <r = 0, (2.2) 

where v denotes the velocity field and cr the Cauchy stress tensor. These statements are the 
conservation of mass and momentum, respectively. For a Newtonian fluid, the extra stress 
tensor T and rate of deformation tensor D are related by 

T = 2r^D, (2.3) 


where rj is a material constant. For an incompressible fluid, the motion of the continuum 
determines the stress tensor up to an arbitrary isotropic tensor and thus a and T are related 
as follows 

«• = -!>I + T, (2.4) 


where p is an arbitrary pressure and I is the identity tensor. 
If we define a stream function ip by 



v = — 


dip 

~di' 


then (2.1) is satisfied identically. Substitution of T from (2.3) into (2.4) and then substitution 
of cr from (2.4) into (2.2) results in the equation 


— Vp + 2pV D = 0. 


(2.5) 


The pressure may be eliminated by taking the curl of (2.5) to give a biharmonic equation 
for the stream function 

Vty = 0. (2.6) 

Consider Stokes flow through the constricted channel defined by |y| = l(x < 0), |y| = 
a(x > 0) and x = 0(a < |y| < 1). The line y = 0 is an axis of symmetry and so only the 
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upper half of the channel needs to be considered. We assume Poiseuille flow at entry and 
exit which, in terms of the stream function, is defined by 

ip(x,y) — > G(y) as x — > — oo, 0 < y < 1, 

ip(x,y) — > G(y/a ) as x — > oo, 0 < y < a 

where G(y) = iy(3 — y 2 ). Along the channel walls we have no-slip boundary conditions. 

The flow region is truncated on entry and exit at distances h and k from the origin, 
respectively. These lengths are chosen to be sufficiently large so that the flow is fully devel- 
oped in the entry and exit sections. In addition, we impose that the normal derivatives of ip 
vanishes at entry and exit. 


3. Domain Decomposition and Spectral Approximation 

The truncated region is divided into three subregions as shown in Figure 1. Within 
each of the subregions, the solution to the biharmonic equation (2.6) is approximated by 
a truncated expansion of Chebyshev polynomials. The solutions are patched by applying 
C 3 continuity conditions in a collocation sense across the subregion interfaces. This is the 
correct order of weak continuity for this problem. With a judicious choice of collocation 
strategy we show that the approximations are pointwise C° and (7 1 continuous across the 
interfaces. 

In region I ip(x,y) is approximated by ip\x,y) where 

Mi Nj 

V ;/ ( X )2/) = G(y) + ^2 ^ amnPm{ X )Qn{y)- (^-^) 

m= 2 n~ 2 


The functions P„{x), Q J n (y ) are suitable linear combinations of shifted Chebyshev polyno- 
mials chosen so that ip\x,y) automatically satisfies the boundary conditions on x = —h(a < 
y < 1) and y = 1 (-h < x < 0). After a little computation we can show that 

P' (x) = T'(x) + *1 r/(x) + /^To'Cz), 2 <m<M u 

where Tj,(x), 0 < m < M/, are the shifted Chebyshev polynomials on [—A, 0] defined by 

T'(x) = T m (^±±), 

and a^, fa are given by 

fa — ( — l) m m 2 , fa — (-l) m (m 2 - 1), 2 <m<M/. 

Similarly, we can show that 

Qi(y) = + t&fHv) + P’Xh), 2 < n < Ni, 

where T^(y), 0 < n < AT/, are the shifted Chebyshev polynomials on [a, 1] defined by 

= T » - 
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and a„, (3^ are given by 


“n = -” 2 , Ph = - 1- 

In region II, the stream function is approximated by where 

Mjj Nn 

V(x,y) = G(y) + E E (3.2) 

m— 2 n= 2 

Since regions I and II share the boundary y = a(—h < x < 0) and we wish the approxi- 
mations to be conforming, we choose Mjj = Mj. Further, since ^>^(x,y) satisfies the same 
boundary conditions along x = -h as ^{x.y) we take P^(x) = P^(x). As for the’ In- 
direction, we define the functions Qlf{y) so as to satisfy the conditions along the axis of 
symmetry. Accordingly Qj/(y) is defined by 

Q!!(y) = n\v) + "(y) + ? "(y), 2 < n < n h , 

where T^(y), 0 < n < Nu, are the shifted Chebyshev polynomials on [0, a] defined by 

«'(») = (^r ) . 

and af/tPn 1 are given by 

& n = ~ 1), = -(-l) n + ^(-l)V(n s - 1), 2 < n < N n . 

In region III, the stream function is approximated by i> HI (x,y) where 


Mm Nui 

«"'(*,») = G(»/o) + E E <*» PL"(x)Qi"(y). (3.3) 

771= 2 71=4 

The functions are suitable linear combinations of shifted Chebyshev poly- 

nomials chosen so that ip II! {x,y) automatically satisfies the boundary conditions on y = 
0(0 < x < A:), y = a(0 < x < A) and x = &(Q < y < a). As before, we can show that 

= #"(*) + «£"*/"(*) + 2 < m < 

where T^ l // (x), 0 < m < Mm, are defined by 


and are g iven b Y 



a" 7 = -m a , PjJ 1 — m 2 — 1. 


/// 

771 


Since we wish to have polynomials of the same order on both sides of the interface x = 0(0 < 
y < a) we choose IV/jj = iV//. The polynomials Qlf\y') are given by 


Qn\y) = fn n (y) + ^ n n Ti n (y) + Pl! l fl n {y) + ti n fl n {y) + S"r 0 7 "(y), 4 < n < A //7 , 
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where fl u (y) = f™(y), 0 < n < N n , and 

= n 1- "’ ' 5 (_1 + (_1) " ) + 5 (_1) " n!( ’ 12 _ *) 1 ' 


S" = 6a"'-i(-l)"nV-D, 

i"' = + i(-i + (-1)’). 

%“ = -1 - a"' - Al" - 7 '". 


4. Collocation Strategy 

The coefficients in the series expansions (3.1) - (3.3) are determined by collocating the 
differential equation at certain points in the domain and the solution in the three subregions 
are patched by imposing the correct order of weak continuity across the subregion interfaces. 
The points at which the extreme values of the Chebyshev polynomials are attained are well 
known to give rise to optimal approximation properties of smooth functions. Therefore, we 
choose the points corresponding to the Chebyshev polynomial of highest degree used in the 
solution representations as our collocation points in both the x- and y-directions. Thus, in 
region I, for example, the collocation points are given by 

j h(xi - 1 ) j (1 - a)yj + 1 + a 

x * ~ 2 ’ Vj ~ 2 

where 

Xi = — cos( 

Vj = -cos(^-), 0 

Boundary conditions 

Due to the choice of modified Chebyshev polynomials as trial functions in the expansions 
(3.1) - (3.3) the boundary conditions are automatically satisfied except along x — 0(a < y < 
1). Along this part of the boundary there are Nj + 1 collocation points. We deduce from 
(3.1) that ip and ip x are polynomials of degree Ni along x = 0(a < y < 1) each depending 
on Ni — 1 degrees of freedom. Therefore, collocation of tp and tp x at Nj — 1 distinct points 
ensures that the boundary conditions along x = 0(a <2/^1) are satisfied identically. We 
collocate these conditions at the points (0, 2//), j = 0, . . . , Ni — 2. 

Interface continuity conditions 

Let us first examine the interface y = a(—h < a < 0) between subregions I and II. We 
impose C 3 continuity of the stream function across this interface, i.e., 

= (*.“)> * = 0 . 1 . 2 , 3 , -h<x< 0 . ( 4 . 1 ) 
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Now ip 1 and ip 11 are polynomials of degree Mj along y = a each possessing Mj — 1 degrees 
of freedom. We may write the condition (4.1) with k — 0 as 


Mil 

E 


n= 2 


E - E *-<?"(«) 


n=2 


n=2 


*>,!.(*) = 0 . 


(4.2) 


We want to collocate at a sufficient number of points to ensure that (4.2) is satisfied identi- 
cally in which case ip is pointwise continuous across the interface between subregions I and 
II. Equation (4.2) is collocated at the points (xf,a),i = 2,3 ,Mj — 2. A further two 
conditions are required to ensure (4.2) is an identity. Therefore, in addition, we collocate 
the following conditions 
Bib 111 

V»"(0,a) = l, —q— (0,a) = 0. (4.3) 

Similarly, continuity of ip y is obtained by collocating (4.1) with k = 1 at the points (x[, a),i = 
2, 3, • • • , Mi — 2, as well as the conditions 


dtp 11 

By 


(0,a) = 0, 


d 2 ip 


n 


dxdy 


(0,a) = 0. 


(4.4) 


This results in pointwise continuity of ip y across the interface y = a(—h < x < 0). 

Now consider the continuity conditions on the second and third derivatives of ip across 
the interface. We collocate each of the conditions (4.2) with k = 2 or k = 3 at the Mj — 3 
points (xj, a),* = 2, 3, — 2. Thus, these derivatives are not pointwise continuous 

across the interface. Moffatt [9] shows that the leading singular term in the expansion of ip 
about the corner (0,a) is 0(r*) where A = 1.5445 and so it would be inconsistent to impose 
continuity of these higher derivatives. 

The same collocation strategy is applied across the interface x = 0(0 < y < a) between 
subregions II and III. As a result the stream function and its normal derivative are pointwise 
continuous across the interface whereas the second and third derivatives of ip are continu- 
ous only at the collocation points (0,yj r ),j = 2, 3 , ••• ,Nji — 2. In addition all tangential 
derivatives of ip and ip y are pointwise continuous across the interface. 


Differential equation 


The biharmonic equation is collocated in each subregion at all points on the collocation 
grid with the exception of those on or one in from each subregion boundary, i.e., 


( x i j 3/j )> ® — 2, • • • , Mk — 2, j — 2, ■ ■ • , Nk — 2, 


for k = I, II, and III. 


When the spectral collocation equations resulting from the boundary conditions, interface 
continuity conditions and differential equations are added together, they yield a total of 
[(Mi — 1 )(Ni — 1) + (Mu — 1 )(Nii — 1) + (Mm — l)(Nm — 3)] linear equations which is equal 
to the number of unknown coefficients a mn , b mn and c^. Therefore, provided the coefficient 
matrix is non-singular, this system of equations possesses a unique solution. 
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The spectral representations we use ((3.1) - (3.3)) are constructed to automatically satisfy 
some of the boundary conditions. It is possible to construct a collocation scheme in which 
the basis functions satisfy none of the boundary conditions but which results in the same 
scheme. This has the advantage of being more flexible but requiring the solution of a slightly 

larger linear system of equations. 

In subregion I, for example, we represent if; in the form 

M N 

«'(*, y) = G(y) + E E -ri(x)^(»). (4.5) 

m=0 n=0 

Satisfaction of the boundary conditions along x = -h(a < y < 1) at the collocation points 
(~h,y]),j = 0, 1, • • • , Ni leads to the following sets of equations: 

Nj ( Mj ) 

E E !?(»/) = 0. 0 <j<N It (4.6) 

71=0 ( 771=0 J 

Nj ( Mj \ 

E E(- i r m2 “™» A'(s') = 0. 0<J< N,. (4.7) 

n=0 ( iTi—O J 

In effect, we are collocating a polynomial of degree Ni Ni + 1 points and therefore since 
we are equating it with the zero polynomial, its coefficients must be zero, i.e., 

M/ 

E(-l) m Omn = 0, 0 <n<N It (4.8) 


Mi 

£) {-l) m m 2 a mn = 0, 0 < n < JV/. 

m— 0 

We may eliminate a 0n and a ln from (4.8) and (4.9) to obtain 

Mi 

a ln = (-l) m ^ 2a mn, 0 <Tl< Nj. 

m=2 


Mj Nj 

'0 / ( X )2/) = G{V) + ^ Zl a mnPm( X )T n {y)> 

m — 2 n=0 


Mj 

a 0n — ZZ (~l) m ( m2 — l) a 7n7ij 
m=2 

Thus, we may write (4.5) in the form 


(4.9) 


(4.10) 


where P^{ x ) * s defined in Section 3. 

The boundary conditions along y = 1 (— h < x < 0) leads to 


Mj 

{ Nj ) 


(4.11) 

E 

\ > 

ft 

II 

J 3 

m= 2 

U=° J 



M Z ( 

' Nj 

1 

(4.12) 

E 


^( x ) = o. 

m—2 \ 

n=0 

J 
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We collocate (4.11) and (4.12) at the collocation points (xj, l),z = 2,3, • • • , M/. As 
before, the coefficients in (4.11) and (4.12) must be zero when collocated at these points, 
i.e., 

Nj 

X] a mn = 0, 2 < m < Mj, (4.13) 

71=0 

Nj 

E n2(1 mn = 0, 2 < m < Mi. (4.14) 

n=0 

Again, we may eliminate a m0 and a ml ,2 < m < Mj, from (4.13) and (4.14) to obtain 

Nj N r 

a m0 2 ' ( n l) a T7inj a mn = ^ ] n Omn- (4.15) 

n=2 n=2 

Substitution of (4.15) into (4.10) results in the approximation (3.1). Thus, we have shown 
the equivalence of the two expansions (3.1 ) and (4.5) in subregion I when an appropriate 
collocation strategy is chosen for the boundary conditions. We may show a similar result in 
the other two subregions. 

5. Method of Solution 

We describe the capacitance matrix technique [1, 15] for efficiently inverting the linear 
system of equations derived in the previous section. It is important to incorporate the 
underlying matrix structure into the solution procedure so that not only do we have an 
efficient method but also one which is able to solve a large problem without running into 
storage problems. 

The spectral domain decomposition method described in this paper gives rise to a block 
matrix whose blocks are either full or zero. The full problem may be solved directly, but this 
process would not be efficient in terms of the number of storage locations and computational 
time. The linear system of equations for the unknown coefficients is written in the partitioned 
form 


o 

O 

0 

H 

1 


x x ' 


’ V 

0 b 2 0 b 4 b & 


x 2 


Q 

<3 

o 

o 

o 

o 


Z3 


0 

D\ D 2 0 D\ D$ 


£4 


0 

o e 2 e 3 f? 4 e 5 




r 


The first three blocks of rows correspond to collocation of the biharmonic equation and 
boundary conditions in subregions I, II, and III, respectively. The last two blocks of rows cor- 
respond to the interface continuity conditions between subregions I and II and subregions II 
and III, respectively. The square matrices A 1} B 2 , C 3 , D 4 , and E 5 are of orders m, n 2 , n 3 , 2n 4 , 
and 2 ri 5 , respectively, where 

m = (Nj - 4 )(Mr - 4) + 2(Ni - 2), 

7*2 = (Nn - 4 )(M n - 4) + 4, 
n 3 = ( Njn — 4 )(Miu — 4), 
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n 4 = 2 (Mj — 4), 
n 5 = 2 (N ni - 4). 

The vectors x* are just a partition of the vector of unknown expansion coefficients in the 
three subregions. 

In the capacitance matrix method, we write (5.1) in the natural component form sug- 


gested by the partitioning 

(5.2) 

A1X1 + A 4 x 4 = p 

B 2 x 2 + B 4 x 4 + B 5 x 5 = q 

(5.3) 

C 3 x 3 + C5X5 = 0 

(5.4) 

D\x 4 + D 2 x 2 + D 4 x 4 + U5X5 = 0 

(5.5) 

E 2 X 2 + E 3 X 3 + E 4 X 4 + E5X5 = T. 

(5.6) 

We write x 4 , x 2 , and x 3 in terms of x 4 and x 5 by premultiplying (5.2), (5.3) and (5.4) by the 


inverses of A\,B 2 , and C 3 , respectively provided they exist, i.e., 

Xi = A* * p — An ^ A 4 x 4y x 2 = B2 q B2 B4X4 B2 B5X5) 

(5 ' ?) 

X3 — — C5X5. 


Eliminating X\, x 2) and X3 from (5.5) and (5.6) we obtain the following system for x 4 and X5. 
(D 4 - D 1 AT 1 A a - D 2 B 2 l B A )x 4 + ( As - D 2 B 2 l B z )xs = -D 1 Az 1 p - D^q, 

(5.8) 

(E 4 - j E 2 j B 2 - 1 S 4 )x 4 + (E 5 - E 2 B 2 1 B a - E 3 C 2 1 C s )x 5 =r- E 2 B 2 1 q. 

The system (5.8) is one of order 2(n 4 + n 5 ) and can be solved efficiently since it is much 
smaller than the original system. The coefficient matrix of this system is known as the 
capacitance matrix. We solve (5.8) simultaneously for x 4 and x$ and then determine Xx,x 2 , 
and X3 from (5.7). Whenever a system of equations needs to be solved in this method, we 
use a Crout factorization technique from the NAG Library [11]. 

The following algorithm describes how efficiencies can be made in the solution procedure. 


Algorithm 

(1) Calculate A^ A a and A^p by solving a system of the form 

A 1 [W 1 \v 1 ] = [A a \p\. (5.9) 

Exploit the fact that A 4 has only n 4 nonzero columns, and therefore we need only solve a 
system with n 4 + 1 right-hand sides. Similarly, we calculate B^ 1 B 4 , B^Bs, B 2 q, and C 3 C$ 
carefully exploiting any zero columns that B A B 5 or C 5 may possess by solving the systems 

B 2 [W 2 \W 3 \v 2 ] = [B 4 \B 5 \q], (5.10) 

G 3 [W 4 ] = [C s ]. (5.11) 


9 


(2) Evaluate the entries in the capacitance matrix: 


D A - D 1 W 1 - D 2 W 2 D 5 - D 2 W 3 

e a - e 2 w 2 e 5 - E 2 W 2 - E 3 W A J ’ 

and the right-hand side: 

— D\V\ — D 2 v 2 
r — E 2 v 2 ’ 

again exploiting any zero columns in W\, W 2 , W 3 and W A . 

(3) Solve the capacitance matrix system (5.8) for x 4 and x 5 . 

(4) Evaluate x 1 ,i 2> and x 3 by substitution of x 4 and x 5 in (5.7). 

The majority of the work in this algorithm lies in Steps (1) and (3). In Step (1) the 
solution of the linear systems (5.9), (5.10), and (5.11) requires Ami-fnJ(n 4 -|-l), &7i2+7i2(27i 4 -f 
715) and knl + 713(715 + 1) operations, respectively. The solution of the capacitance matrix 
system requires 8k(nl + 715) + 4(71^ + n 5) operations. The solution of the full system (5.1) 
would require kN 3 + N 2 operations where N = ni +n 2 + 7x3 + 2(714-1-715). Thus, the use of the 
capacitance matrix technique has effected savings of O^i^+r^+n^nx+raa) +713(711+713)). 

6. Numerical Results 

Numerical experiments are performed for Stokes flow through a 4:1 contraction geometry, 
i.e., o = l. We examine the convergence of the approximations as the degree of the trial 
functions is increased and also the performance of the capacitance matrix method for solving 
the spectral collocation equations. As in [6, 15], we choose truncation parameters h = 1.5 
and k = 0.5. 

To verify the convergence of the approximations (3.1) - (3.3), we give contour plots of 
the stream function for different numbers of degrees of freedom. In all of these plots, we 
have Mi = Mu = Mm = 14. In the y-direction, we choose AT/ = Nu — 12, AT/ — AT// — 16 
and Ni = Nu = 20 in Figures 3, 4, and 5, respectively, with AT/// = 8 in all cases. The 
smoothness of the contours improves as the number of degrees of freedom is increased. Notice 
also the continuity of the stream function across the subregion interfaces which is a result of 
the conforming discretization that we use. 

The contours appear a little ragged, especially those defining the vortex in the salient 
corner. This is due to the contouring routine. Although we obtain a global approximation, 
the NAG routine for contouring a continuous function fails because of sharp charges in the 
gradient of the solution in the recirculation region. To overcome this, a uniform mesh is 
placed over the domain and a NAG routine interpolating the values of the solution at these 
points is used. 

The stream function plot in Figure 5 obtained with a total of 571 degrees of freedom is 
in qualitative agreement with those in [8] which are obtained using a finite element method 
with as many as 1326 degrees of freedom. A significant improvement is also observed over the 
authors’ previous work [5, 15] in the continuity achieved across the subregion interfaces which 
indicates the advantage of using a conforming spectral collocation strategy. The contour plot 


E 

} 
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of the vorticity is given in Figure 6. This is obtained by differentiating the stream function 
shown in Figure 5. 

In Table I, we compare the direct inversion of (5.1) with the capacitance matrix method 
in terms of CPU (Amdahl 5890 - 300) time and required storage in Mbytes for different 
numbers of degrees of freedom. In this comparison, we choose M/ = Mjj = Mm = Nj = 
Nu = N m. From the table, we note that the capacitance matrix method becomes more 
efficient as the number of degrees of freedom increases in terms of both storage requirements 
and computational effort. Alternatively, with a given storage limit, one may solve, using 
the capacitance matrix method, a problem using many more degrees of freedom than by the 
original method which inverts directly the linear system (5.1). 

7. Conclusions 

A spectral domain decomposition method is described for solving Stokes flow in a channel 
contraction using the stream function formulation. Spectral approximations are constructed 
and a collocation strategy devised so that the subregions are conforming. The resulting ap- 
proximations and their normal derivatives are pointwise continuous across subregion inter- 
faces. The trial functions in these representations are chosen to satisfy some of the boundary 
conditions. An alternative collocation method in which the trial functions do not satisfy any 
of the boundary conditions is shown to be equivalent to the former strategy. 

Efficient direct methods based on the capacitance matrix technique are used to solve the 
resulting system of linear equations for the expansion coefficients in the three subregions. 
Economies are made in terms of storage and computational effort over the original method. 
This solution procedure can be applied to any rectangularly decomposable domain. 

In previous work [6], the use of nonconforming subregions in the solution of the Navier- 
Stokes equations led to a breakdown in convergence at Reynolds numbers above 200. We aim 
to use the new conforming subregions to extend the range of Reynolds numbers for which 
the method converges and this will be reported in a future paper. 

Table 1. CPU (Amdahl 5890-300) times in seconds and storage in Mbytes required for 
solving the 4:1 contraction problem for different numbers of degrees of freedom. 


Mi 

Degrees of Freedom 

Original Method 

Capacitance Matrix 



time 

mbytes 

time 

mbytes 

9 

176 

1.35 

0.97 

1.44 

0.94 

11 

280 

2.57 

1.35 

1.99 

1.19 

13 

408 

5.89 

2.06 

3.40 

1.54 
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Figure 1. Three subregion decomposition of the flow region. 



Figure 2. Two subregion decomposition of the flow region. 
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Figure 6. Vorticity contours of the 4:1 problem with 571 degrees of freedom. 
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